The  Statistics  of  GPS 


Demetrios  Matsakis 

U.S.  Naval  Observatory,  3450  Massachusetts  Avenue  NW,  Washington,  DC.,  USA,  20392 

ABSTRACT 

The  Global  Positioning  System  (GPS)  is  an  extremely  effective  satellite-based  system  that  broadcasts  sufficient 
information  for  a  user  to  determine  time  and  position  from  any  location  on  or  near  the  Earth.  The  fundamental  GPS 
measurement  is  the  corrected  time  of  the  satellite  clock  relative  to  the  receiver  clock.  This  paper  uses  publicly  available 
information  to  present  a  statistical  analysis  of  the  underlying  timescale  and  clock  performance,  which  can  be  largely 
presented  without  recourse  to  the  many  significant  and  interesting  scientific  corrections  and  parameterized  models  that 
could  or  must  be  applied  to  the  data 
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1.  INTRODUCTION 

The  international  standard  for  time  is  Coordinated  Universal  Time  (UTC),  which  is  published  monthly  in  the  Circular  T 
by  the  International  Bureau  of  Weights  and  Measures  (BIPM)1.  UTC  is  realized  in  real-time  by  the  master  clocks  of 
each  participating  laboratory  k,  and  designated  UTC(k).  The  U.S.  Naval  Observatory’s  (USNO)  realization  of  time  is 
designated  UTC(USNO).  GPS  broadcasts  corrections  so  that  all  of  its  clocks  can  provide  a  predictive  realization  of 
UTC(USNO). 

The  GPS  corrections  are  obtained  with  the  use  of  data  measured  from  what  are  now  1 1  monitor  sites  distributed  fairly 
evenly  in  longitude,  and  located  in  the  non-arctic  latitudes.  Each  monitor  site  records  the  difference  between  the  time 
broadcast  on  each  satellite’s  signals  with  the  time  of  its  local  clock.  These  differences  are  corrected  for  many  effects 
beyond  the  scope  of  this  paper  (such  as  leap  seconds  which  are  intermittently  inserted  into  UTC  to  account  for  the 
rotation  of  the  Earth2),  and  used  by  the  GPS  Master  Control  Station  as  input  to  a  Kalman  filter3  which  models  and 
predicts  satellite  clock  performance  in  discrete  15-minute  steps4.  Measurements  from  the  USNO  are  also  applied  to 
steer  an  implicit  ensemble  mean  to  approximate  a  prediction  of  UTC(USNO) 5. 

The  clock  model  applied  by  the  GPS  Master  Control  Station  Kalman  Filter  (MCSKF)  considers  the  phase  (time  can  be 
considered  to  be  the  phase  of  a  sinusoidal  oscillation),  frequency,  and  frequency  derivative  (frequency  drift)  to  be 
parameters,  and  therefore  this  paper  will  first  describe  the  Allan  and  Hadamard  statistics6,  which  are  appropriate  under 
the  circumstances.  The  paper  will  then  describe  general  considerations  for  a  timescale,  and  provide  details  of  how  the 
GPS  composite  clock,  whose  realization  is  termed  GPS  Time,  operates.  There  are  many  ways  the  GPS  composite  clock 
could  be  steered  to  UTC(USNO),  and  this  paper  will  show  one  way  to  change  the  current  algorithm. 

2.  THE  ALLAN  AND  HADAMARD  DEVIATIONS  AND  VARIANCES 

Statistical  measures  used  in  precise  timekeeping  are  nonstandard  in  other  fields.  Fourier  techniques  are  rarely  used 
because  periodic  terms,  such  as  diurnal  cycles,  are  frequently  minimized  through  environmental  isolation.  Since  no 
clock  property  is  statistically  stationary  over  long  averaging  or  sampling  times,  statistical  measures  that  depend  upon 
differences  over  finite  intervals  of  length  x  are  preferred.  The  Allan  deviation  is  the  most  commonly  used  statistical 
measure,  which  we  define  in  this  section  along  with  the  Hadamard  deviation. 

To  define  these  statistics,  we  denote  the  time  reading  of  a  clock  by  xk(f)  where  f  is  the  time  the  clock  is  measured.  The 
f  are  assumed  to  be  measured  at  equal  discrete  intervals  x.  The  average  frequency  y*  of  a  clock  between  f  and  ti+i  is 
given  by: 


yi=(xi+1-xi)/x  (1) 
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The  Allan  deviation,  denoted  ay(x),  is  the  square  root  of  the  Allan  variance,  a2y(x),  which  is  a  scaled  mean  square  of  the 
second  differences  of  the  clock  time  readings,  or  the  first  difference  of  the  frequencies: 

a2y(x)  =<(xi+2-  2  xi+i  +  Xi)2>  1(2  t2)  (2) 

Note  that  this  statistic  is  independent  of  the  clock’s  initial  time  and  frequency.  The  division  by  two  normalizes  the 
statistic  so  the  Allan  deviation  would  equal  the  RMS  of  the  clock  that  is  characterized  by  white  Gaussian  noise. 

The  Hadamard  deviation,  denoted  Hc>y(T),  is  the  square  root  of  the  Hadamard  variance,  H^2yC0?  which  is  a  scaled  mean 
square  of  the  third  difference  of  the  time,  or  the  second  difference  of  the  frequencies: 

H<J2y(x)  =<(xi+3  -  3  xi+2  +  3xi+1  -  Xi)2>  /(6  T2)  (3) 


Note  that  this  statistic  is  independent  of  the  clock’s  average  time,  frequency,  and  frequency  drift.  It  is  sensitive  to  any 
change  in  the  frequency  drift,  or  to  any  change  in  the  time  or  frequency  which  cannot  be  characterized  through  an 
unchanging  frequency  offset  or  frequency  drift.  The  reader  will  have  already  surmised  that  the  motivation  for  this 
statistic  was  to  characterize  clocks  with  unpredictable  initial  and  average  values  of  the  time,  frequency,  or  frequency 
drift.  These  are  the  rubidium-based  atomic  clocks  which  are  currently  used  in  the  GPS  satellites7. 

3.  CONSIDERATIONS  IN  GPS  TIMESCALE  FORMULATION 

A  timescale  can  be  described  as  a  weighted  average  of  corrected  clocks.  Timescales  can  be  useful,  if  not  required, 
because  there  is  no  way  to  measure  absolute  time8,  there  is  no  perfect  clock,  and  because  only  clock  differences  are 
measurable.  For  many  applications,  it  is  essential  that  the  timescale  be  smoothly  varying,  and  robust  with  respect  to  the 
addition,  subtraction,  or  reweighing  of  individual  clocks.  In  most  cases,  a  real-time  timescale  is  defined  through 
extrapolation  of  a  previously-determined  timescale;  the  current  value  of  the  timescale  is  set  equal  the  weighted  average 
of  the  difference  between  the  measured  and  predicted  clock  values. 

Although  the  averaging  process  can  easily  create  a  timescale  which  at  least  equals  the  performance  of  an  individual 
clock  with  regards  to  any  statistic  of  interest,  it  is  not  possible  to  create  a  timescale  that  is  simultaneously  optimal  with 
respect  to  every  statistic.  For  example,  the  goal  of  the  GPS  system’s  composite  clock  is  to  create  a  timescale  that  is 
optimal  for  navigational  solutions,  but  not  necessarily  to  optimize  the  extraction  of  UTC(USNO).  Therefore,  a 
correction  is  broadcast  with  the  message,  so  that  UTC(USNO)  can  be  extracted  via  a  deterministic  correction  to  GPS 
Time  (Figure  l)5.  The  compromises  involving  in  steering  GPS  time  will  be  discussed  in  part  V. 
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Figure  1.  Upper  plot  is  the  daily  average  of  UTC(USNO)-GPS  Time,  as  measured  at  the  USNO.  Lower  plot  is 
the  difference  between  UTC(USNO)  and  GPSfs  delivered  prediction  of  UTC(USNO).  Data  are  in  nanoseconds, 

and  MJD  is  the  number  of  days  since  November  18, 1858. 


The  question  of  how  to  weight  clocks  in  a  timescale  would  be  easily  handled  in  the  case  of  stationary  Gaussian 
statistics,  but  in  practice  is  complicated  by  the  non-Gaussian,  non-stationary,  and  not  always  measurable  statistical 
nature  of  clocks9.  Even  a  phenomenological  approach  that  characterizes  a  clock’s  statistics  against  the  timescale 
requires  special  care  because  a  clock’s  contribution  to  the  timescale  has  the  effect  of  lowering  the  difference  between 
the  clock  and  the  timescale.  The  stronger  correlation  would  then  increase  that  clock’s  weight  in  the  timescale.  Because 
every  clock  is  100%  correlated  with  its  contribution  to  the  timescale,  this  would  then  lower  the  difference  still  more  as 
the  timescale  is  iterated  or  propagated.  Ultimately,  one  clock  would  carry  all  the  weight  in  the  timescale.  When  one 
clock  carries  all  the  weight  by  design,  the  timescale  is  sometimes  said  to  follow  the  “master  clock  approach”.  For  a 
short  time  during  the  GPS  development  phase,  the  MCSKF  selected  a  single  ground  station  clock  for  this  purpose.  No 
parameters  were  used  by  the  MCSKF  to  describe  this  one  clock’s  state,  and  all  other  clock  parameters  effectively 
described  the  difference  between  the  clock  in  question  and  the  GPS  master  clock.  The  problems  with  that  approach 
became  very  evident  when  that  master  clock  failed,  and  subsequently  the  MCSKF  adopted  the  composite  clock 
formulation  that  is  described  next. 

To  describe  each  system  clock,  three  parameters  per  clock  are  employed,  which  estimate  its  time,  frequency,  and 
frequency  drift.  The  composite  clock  (cc)  is  an  implicit  weighted  average  of  the  corrected  clocks  (the  average 
difference  between  the  clock  readings  and  their  modeled  readings).  The  cc  is  realized  by  the  corrected  value  of  every 
clock  in  the  GPS  system,  but  it  is  never  actually  computed  as  a  numerical  quantity  because  doing  so  is  not  needed  to 
propagate  the  clock  states  in  the  Kalman  formalism.  The  MCSKF  propagates  each  clock  state  forward  based  only  upon 
the  measured  differences  of  the  satellite  and  monitor  station  clocks.  As  an  additional  step,  all  GPS  satellite  and  ground 
clock  frequency  drift  parameters  are  incremented  equally  on  the  basis  of  independent  ground  measurements  of 
UTC(USNO)-GPS  by  the  USNO.  As  shall  be  explained  below,  this  steering  keeps  the  cc  stable  and  in  line  with 
UTC(USNO)  on  timescales  of  several  days  and  more.  Therefore  estimates  of  the  process  noise  that  are  based  upon 
long-term  averaging  of  observations  are  relative  to  UTC(USNO)  and  thus  insensitive  to  the  weights  (inverse  of  the 
squared  uncertainties)  applied  by  the  GPS  MCSKF.  New  clocks  are  introduced  to  the  constellation  with  zero  weight, 
and  their  state  parameters  are  set  to  correct  them  to  the  cc.  Once  weighted,  the  weight  of  any  new  clock  would  increase 
as  more  measurements  are  applied  until  it  reaches  the  natural  level  set  by  the  system  process  noises,  which  are  not 
frequently  changed.  Individual  observations  (of  clock  differences)  perturb  the  two  relevant  clocks  in  opposite 
directions  and  in  accordance  with  the  clock  weights,  therefore  GPS  Time  varies  with  the  weighted  average  of  its  clock 
readings  and  the  MCSKF  solutions  are  relevant  only  insofar  as  they  affect  the  clock  weights4 10 . 

Since  UTC  is  external  to  the  GPS  clocks,  steering  GPS  Time  to  UTC(USNO)  is  essential.  This  is  achieved  through 
acceleration;  the  frequency  drift  parameter  of  every  clock  is  adjusted  by  either  0  or  ±1.0-10"19/s.  A  one-sentence 
summary  of  the  algorithm  is  that  the  acceleration  is  set  opposite  to  the  slope  in  GPS-UTC(USNO),  except  that  if  doing 
so  constantly  thereafter  would  result  in  the  extrapolation  of  the  accelerated  GPS-UTC(USNO)  never  becoming  zero,  in 
which  case  the  acceleration  has  the  same  sign  as  the  slope11.  Figure  2  illustrates  the  process. 
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Figure  2.  Pictorial  display  of  GPS  steering  under  three  sample  circumstances.  Day  0  is  the  time  of  the  steer 

computation,  and  day  -1  is  the  previous  day. 


This  steering  in  acceleration  is  often  termed  bang-bang,  because  it  is  quantized  as  on  or  off.  The  acceleration-based 
steering  creates  a  tendency  for  GPS  timescale  to  develop  the  oscillations  evident  in  the  upper  half  of  Figure  1. 
Historically,  the  oscillations  in  GPS  Time  can  be  described  as  quasi-periodic  with  a  typical  duration  of  20  days.  One 
motivation  for  selecting  this  form  of  steering  would  be  to  attempt  to  minimize  the  subdaily  variations  in  GPS  Time, 
whose  realization  in  the  corrected  satellite  clocks  is  dependent  upon  uploads  that  are  frequently  delivered  no  more 
rapidly  than  once  a  day  per  satellite.  A  situation  frequently  exists  wherein  the  receiver  is  simultaneously  observing 
satellites  whose  broadcast  time  corrections,  including  the  applied  steering,  are  up  to  one  day  out-of-date  compared  to  the 
current  state  values  of  the  MCSKF.  Since  GPS  Time  is  designed  for  position  solutions  that  require  all  satellites  to  be 
on  a  common  time,  but  not  necessarily  the  correct  time,  the  important  thing  is  to  optimize  stability  on  the  shortest 
averaging  times.  The  accompanying  long-term  oscillations  were  deemed  acceptable  because  corrections  to  recover 
UTC(USNO)  are  broadcast  in  the  satellite  navigation  message  so  as  to  generate  the  data  whose  daily  averages  are 
shown  in  the  lower  plot  of  Figure  1. 


IV.  OBSERVED  CLOCK  AND  GPS  TIME  STATISTICS 


The  International  GNSS  Service  (IGS)12  uses  advanced  carrier-phase  techniques  to  compute  the  GPS  satellite  orbits, 
clocks,  ionosphere  corrections,  and  other  parameters.  As  an  associate  analysis  center  of  the  IGS,  the  USNO  provides 
clock  solutions,  which  can  be  used  to  describe  the  difference  between  GPS  satellite  clocks  and  UTC(USNO).  Figure  3 
shows  the  observed  Hadamard  deviation  of  the  GPS  satellite  clocks,  using  USNO-reduced  data  from  a  recent  20-day 
period,  as  reported  to  the  IGS.  The  receiver  was  referenced  to  UTC(USNO),  whose  time  reference  is  considerably  more 
stable  than  any  of  the  GPS  satellite  clocks. 


Hadamard  Deviation,  h<t  (x)  reference  rcvr=  USNO 


Figure  3.  Statistical  characterization  of  GPS  satellite  clocks  as  measured  by  the  receiver  USNO,  which  is 
referenced  to  UTC(USNO),  MJD  54140-54160. 

Most  of  the  observed  GPS  satellite  stabilities  are  identical  to  within  a  factor  of  2,  although  a  few  satellites  are  more  than 
an  order  of  magnitude  worse  in  stability.  An  algorithm  which  created  GPS  Time  by  averaging  all  satellites  equally 
would  generate  results  corrupted  by  the  contribution  of  the  few  noisy  satellites,  such  as  SV15  in  the  figure.  In  the  event 
of  a  satellite  clock  failure,  the  GPS  timescale,  as  seen  by  the  user,  would  show  such  effects  before  a  bad-data  message 
can  be  transmitted  to  the  satellite.  Here  we  discuss  only  timescales  which  apply  weights  according  to  some  statistical 
measure  of  the  clock  stabilities,  and  ignore  the  data-integrity  algorithms  applied  by  the  MCSKF.  One  could  use  the 
Hadamard  deviations  shown  in  Figure  3  to  compute  the  expected  Hadamard  deviations  of  a  timescale  that  weighs  each 
clock  by  its  inverse  Hadamard  variance.  Such  a  timescale  could  base  its  weights  on  the  statistics  at  any  one  of  the 
intervals  (x  ),  and  therefore  be  optimal  for  that  interval  but  not  necessarily  for  other  intervals,  The  computation  of  the 
expected  timescale  stability  is  a  straightforward  procedure  analogous  to  more  standard  mathematics  which  would  apply 
to  a  least  squares  evaluation  involving  the  variances  of  a  normal  distribution.  Figure  4  shows  the  expected  timescale 
stabilities,  and  also  shows  the  observed  Hadamard  deviation  of  GPS  Time  using  the  operational  single-frequency  SPS 
receiver  at  the  USNO. 


GPS  Time  Stability 


Figure  4.  Stability  of  GPS  Time  (least  negative  curve)  plotted  against  the  theoretical  limitations  (lower  plots). 

The  discrepancies  are  explained  in  the  text. 

Because  the  individual  curves  of  the  clocks  plotted  in  figure  3  rarely  cross,  the  differences  between  the  optimization 
schemes  are  barely  noticeable.  However,  the  measured  value  of  GPS  Time  instability  is  considerably  larger  than  would 
naively  be  expected  by  the  clock  stabilities,  particularly  over  the  smaller  averaging  times.  This  is  due  to  the  many 
neglected  sources  of  measurement  noise  in  the  single- frequency  GPS  receiver,  which  does  not  use  carrier-phase 
observations  and  relies  upon  ionosphere,  orbit,  and  clock  model  parameters  that  are  considerably  less  accurate  when 
broadcast  than  the  model  parameters  that  are  current  inside  the  MCSKF.  However,  it  is  evident  that  the  discrepancies 
become  less  as  their  averaging  time  approaches  and  exceeds  a  day,  and  this  is  strong  evidence  that  the  GPS  composite 
clock’s  accuracy  is  not  far  from  its  theoretical  limits  at  large  averaging  times. 

V.  STEERING  GPS  TIME 

Most  precise  timing  applications  involve  steering  clocks  in  frequency,  rather  than  the  frequency  drift  (acceleration) 
steering  used  by  GPS.  In  a  typical  situation,  a  frequency  steer  is  applied  that  is  described  by  the  dot  product  of  a  2  (or 
more)  component  gain  function  vector  with  a  clock  state  offset  vector  that  gives  the  time  and  the  frequency  offset  (and 
perhaps  higher  order  parameters)  of  the  clock  being  steered.  It  is  possible  to  compute  an  optimal  gain  function  for  any 
given  “cost  function”,  which  assigns  a  relative  importance  of  stability  in  time,  stability  in  frequency,  and  the  amount  of 
control  effort,  by  using  a  Linear  Quadratic  Gaussian  (LQG)  formalism13"14.  It  has  been  shown  that  LQG  techniques  can 
be  used  to  generate  a  proportional  steering  method  that  would  enable  closer  steering  of  GPS  Time  to  UTC(USNO),  for 
both  the  existing  and  expected  future  GPS  constellation  parameters15.  Figure  5  and  6,  taken  from  that  work,  show  how 
the  time  determinations  could  be  improved  by  proportional  steering  in  a  future  GPS  constellation. 
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Figure  5.  Simulated  performance  of  a  future  GPS  constellation  using  the  current  GPS  acceleration-based 
steering  algorithm 
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Figure  6.  Simulated  performance  of  cesium-based  GPS  constellation  using  a  proportional  steering  algorithm 
selected  to  minimize  frequency  instabilities  as  opposed  to  time  instabilities. 


Figure  7  uses  the  Allan  deviation  to  quantify  the  improvement15.  This  would  be  the  appropriate  statistic  for  a 
constellation  using  cesium  atomic  frequency  standards,  or  any  clocks  that  did  not  require  estimation  of  the  frequency 
drift. 


Tau(sec) 

Figure  7.  Allan  deviations  associated  with  simulations  of  previous  two  figures. 


SUMMARY 


The  statistical  properties  of  the  GPS  system  can  be  characterized  by  standard  precise  timekeeping  measures,  such  as  the 
Allan  and  Hadamard  deviations.  Future  work  involving  clock  steering  could  prove  fruitful,  particularly  as  the  GPS 
system  and  performance  improve,  and  to  improve  interoperability  with  other  GNSS  systems. 
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